function Qdraws = drawQs_ar(restr,phi,opt)
% Function draws Qs from the space of orthonormal matrices 
% satisfying sign, elasticity and narrative restrictions.
% Function assumes identified set is nonempty (if empty, function will
% never terminate).
% Inputs:
% - restr: structure containing information about restrictions
% - phi: structure containing reduced-form VAR parameters
% - opt: structure containing model information and options

sphi = restr.sphi;
hdSignRestr = restr.hdSignRestr;
Sigmatrinv = phi.Sigmatrinv;
vma = phi.vma;
U = Sigmatrinv*restr.U; % Multiply by Sigmatrinv to avoid doing this repeatedly
Sigmatrinvp = Sigmatrinv'; % Transpose to avoid doing this repeatedly

n = size(Sigmatrinv,1); % Number of variables in VAR
mHD = size(hdSignRestr,1); % No. of restrictions on historical decomp.

%% Draw Qs satisfying restrictions.
Qdraws = zeros([n,n,opt.qDraws]);

parfor kk = 1:opt.qDraws % Replace with for if no Parallel Computing Toolbox

    %[kk opt.qDraws]

flag = 0;

while flag == 0

    % Draw Q from space of orthonormal matrices.
    z = randn(n);
    [Qtilde,~] = qr(z);    

    % Normalise diagonal elements of A0 to be positive. Note that Matlab 
    % is implicitly expanding arrays to be compatible with elementwise 
    % array operations.
    Q = ((sign(diag(Sigmatrinvp*Qtilde))').*Qtilde)./vecnorm(Qtilde);

    % Check whether proposed draw satisfies sign, elasticity and shock-sign
    % restrictions.
    sSat = all(sphi*Q(:) >= 0);
       
    if sSat == 0  
        % If restrictions not satisfied return to beginning of loop.
        continue;          
    end

    % Check whether proposed draw additionally satisfies restrictions on 
    % historical decomposition.
    hdCheck = zeros(mHD,1); 

    for ii = 1:mHD % For each restriction  
    
        % Pre-compute q_j*q_j' for j=1,...,n.
        QQ = pagemtimes(reshape(Q,[n,1,n]),reshape(Q,[1,n,n]));

        % Extract row of hdSignRestr corresponding to ith restriction.
        hdRestr = hdSignRestr(ii,:);

        % Compute contribution of each shock.
        HH = zeros(1,1,n);    

        for hh = 1:hdSignRestr(ii,4)+1 % Sum contribution over horizons

            HH = HH + pagemtimes(pagemtimes(vma(hdRestr(1),:,hh),QQ),...
                U(:,hdRestr(3)+hh-1));

        end

        if hdRestr(5) == 1 && hdRestr(6) == 1

            % Type A - most important contributor
            hdCheck(ii) = abs(HH(hdRestr(2))) == max(abs(HH));

        elseif hdRestr(5) == 1 && hdRestr(6) == -1

            % Type A - least important contributor
            hdCheck(ii) = abs(HH(hdRestr(2))) == min(abs(HH));

        elseif hdRestr(5) == 2 && hdRestr(6) == 1

            % Type B - overwhelming contributor
            hdCheck(ii) = abs(HH(hdRestr(2))) - ...
                sum(abs(HH((1:n) ~= hdRestr(2)))) >= 0;                

        elseif hdRestr(5) == 2 && hdRestr(6) == -1

            % Type B - negligible contributor
            hdCheck(ii) = abs(HH(hdRestr(2))) - ...
                sum(abs(HH((1:n) ~= hdRestr(2)))) <= 0;

        end
        
        if hdCheck(ii) == 0              
            % If ith restriction on historical decomposition not 
            % satisfied, exit loop (saves computing contributions 
            % relating to subsequent restrictions).
            continue;
        end

    end

    hdSat = all(hdCheck);

    if sSat*hdSat == 1

        flag = 1;
        Qdraws(:,:,kk) = Q;

    end

end

end